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1. Introduction 

The first solution of a finite-dimensional lattice model that exhibits a phase 
transition was the calculation by Lars Onsager in 1944 of the partition function of 
the zero-field square-lattice Ising model. [I] In 1949 Bruria Kaufman derived this 
result using spinor or free-fermion operators. [2] 

We discuss the history of the derivation of the spontaneous magnetization Ai 
of the Ising model in section 4: it depended on first using the free fermion structure 
to write A4, for a finite lattice, in terms of a determinant, then taking the limit 
when the lattice and the dimensions of the determinant become infinite. [3j HJ [5] 

After the Ising model, and the related dimer problems, the next models to be 
solved were the six-vertex and eight-vertex models. We discuss these in sections 5 
and 6, in particular we outline the Bethe ansatz calculation for the general six-vertex 
model in a field. Then many other solvable models were found, both two and three- 
dimensional. [6] -[15] These were more complicated than the Ising model, lacking its 
simple free-fermion structure and explicit determinantal expressions for correlations 
such as the magnetization. (One exception is the six-vertex model with particular 
boundary conditions. [To]) 

One of the most challenging has been the two-dimensional iV-state solvable 
chiral Potts model. [17] It lacks the "rapidity-difference" property, so there is no 
explicit single-valued parametrization of the Boltzmann weights of the model, and 
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the simple inversion relation and corner transfer matrix tricks to calculate the free 
energy and magnetization fail. 

However, there is a special "superintegrable" case [18] of the chiral Potts model 
whose magnetization is that of the general model and which has simple properties 
similar to those of the Ising model. Indeed, when N = 2 it is the Ising model. 

In particular it has recently been shown [T9]- [23] that its magnetization can be 
expressed as an m by m determinant D, even for a finite lattice. 

So the wheel has come full circle and we are back to determinantal expressions. 
It still remains to calculate D in the large-lattice limit, when m becomes infinite. 
For the superintegrable chiral Potts model, the determinant is not Toeplitz, but it is 
Cauchy-like, so D can be evaluated explicitly and the magnetization obtained. [22] 

2. Partition functions and transfer matrices 

We draw the square lattice diagonally, as in Figure [U with M rows of sites 
and L sites per row. On each site i we place a "spin" <7j, which takes some set of 
discrete values. For the Ising model, Oi = +1 or — 1; for a general iV-state model, 
o"j = 0, 1, . . . , N — 1. Adjacent spins i,j interact with an energy e(<7j, Oj). Spins on 
the top and bottom rows are fixed to the value (for the Ising model, to the value 
+ 1), and cylindrical boundary conditions are imposed, so that we identify the spins 
in the extreme left column with those in the extreme right. 

Let 

W{ffi,<Tj) = e - e{a ^ )/kBT 

be the Boltzmann weight of the edge {ks is Boltzmann's constant and T the 
temperature). Then the partition function is 

z = £ Uwfatrj) , (2.i) 

where a is the set of all spins, and the sum is over all their permitted values. 
We are particularly interested in calculating the partition function per site 

k = Z^ LM , (2.2) 

the dimensionless free energy f = — log k, and averages such as the Ising model 
magnetization 

M = Z- x Y,°i • ( 2 - 3 ) 

(hi) 

We expect k and M. to tend to limits when L, M — > oo and the central spin o\ 
becomes deep within the lattice. 

Of course, such calculations are not easy: for the Ising model on a ten by ten 
lattice there are 2 100 ~ 10 30 terms in the sum in f)2.ip . One starts by defining the 
row-to-row transfer matrix T. Let s = {a±, a 2 , ■ ■ ■ be the spins in one row of 
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Figure 1. The square lattice turned through 45°. 



the lattice, and s' = {a[, o~' 2 , . . . , cr' L } the spins in the row above. Define the N L by 
N L matrix T with entries 



t.s = n%^)%/i-ii 



(2.4) 



taking W to be the weight on SW-NE edges and W to be the (possibly different) 
weight on SE-NW edges. Then the partition function is 

Z = v)T M u , (2.5) 

where u is the iV L -dimensional vectors with entries 

u s = S(ax,0)---5(a L ,0) . (2.6) 

For M large, it follows that 

Z ~ (A max ) M , (2.7) 

so we have reduced the calculation of k to the calculation of the maximum eigenvalue 
of an A rL -dimensional matrix. For the above-mentioned N = 2, L = 10 case this 
dimension is 2 10 = 1024, which is a huge improvement on 10 30 ! 
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3. Ising model 

Even so, one still wants to take the limit L — > oo, so T is going to become 
infinite-dimensional. For the Ising model, a, t = ±1 and 

Wfaaj) = e K ^ , Wfaaj) = . (3.1) 

Onsager considered this model, with the more usual 90° orientation of the lattice and 
toroidal (cyclic) boundary conditions. Then T is the product of two matrices, one of 
which adds the horizontal edges within a row, while the other adds the vertical edges 
between rows. He showed that these two matrices generate a finite-dimensional Lie 
algebra, now known as the Onsager algebra. [H eqns. 60 & 61]. This enabled him to 
calculate all the eigenvalues of T, and hence Z and k. 

The matrix T commutes with the operator R that negates all spins: 

L 

Rs,s> = USfa-tf) , R 2 = l (3.2) 

i=i 

so its eigenvectors v either lie in the sub-space where Rv = v or in the sub-space 
Rv = —v. In the former case they are contained in the eigenvalues of T + , in the 
latter in T_, where 

T± = t 1 ®t 2 ®---®t L (3.3) 

and the ti are two- by-two matrices (the are different for T + and T_). So in this 
sense T is a direct product of L two-by-two matrices. 

Bruria Kaufman later gave a simpler derivation of this result, using anti- 
commuting spinor (free-fermion) operators, i.e. a Clifford algebra. [2] 



4. Ising model Magnetization 

4-1. Kaufman and Onsager 

In 1949, at a conference in Florence, Italy, [3] Onsager referred to the 

magnetization of the Ising model and announced that "B. Kaufman and I have 
recently solved" this problem. He gave the result as 

M = (1-k 2 ) 1 / 8 (4.1) 

where 

k = l/(sinh IK sinh 2K) (4.2) 

and the result is true for < k < 1. For k > 1 the magnetization vanishes, i.e. 
Ai = 0, so the graph of Ai is as in Figure [2j 

Onsager and Kaufman did not publish their derivation, but we have convincing 
evidence of their method. In 1949 they published a paper[23] entitled "Crystal 
statistics. Ill" where they use the free-fermion spinor operators to calculate the 
correlation (si^s^i+fc) for i — 1, 2 and any integer k [not the k of (14.21) ) in terms of 
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Figure 2. A4 as a function of temperature T . 

a /c-dimensional Toeplitz determinant A*.. The magnetization can be obtained from 
these expressions by taking the limit k — > oo, when (s lt is ijl+k ) — > M 2 . 

So they needed to calculate the large-size limit of a Toeplitz determinant. In 
a letter dated April 12, 1950, Onsager wrote to Kaufman giving a general formula 
for this limit Aoo.[25] On page 3 he says "we get the degree of order from C.S.III 
without much trouble. It equals (1 — k 2 ) 1 ^ 8 as before". This letter is also in the 
Onsager archive in Trondheim, at pages 21 - 24 of 
http: / / www.ntnu.no /ub/ spesialsamlingene/tekark/tek5 / research/009_0097.html 
In the same archive, on pages 32, 33 of 

http: / / www.ntnu.no/ub / spesialsamlingenc /tckark / tck5 / rcscarch/009_0096.html 
there is a letter from Kaufman to Onsager thanking him for his letter and saying 
she has worked out a way of using long-range order along a row, rather than a 
diagonal, and applying his procedure to obtain A^. She goes on to say that the 
mathematician Kakutani had written to her saying that he had spoken to Onsager 
about this, and was very interested. 

Onsager recounts something of what happened in 1971. In [26J he says that 
Kaufman derived the correlations in the form of recurrent (Toeplitz) determinants, 
and they were particularly simple along a diagonal (i.e. a row or column of Figure 
[1]). He then discusses the problem of calculating the infinite-size limit of these 
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determinants and indicates how he first did this for the particular problem by 
determining the eigenvalues and taking their product. This was the basis for his 
announcement of the result in Florence. 

He then looked for a more general formula for Toeplitz determinants with 
arbitrary generating functions. He found one for a large class of rational generating 
functions. 

From this, he says "the general result stared me in the face" . This is the formula 
(7) of [25], which is exact for for his rational generating functions when the dimension 
k of the determinant is finite but sufficiently large. Since the elements of are the 
central 2k — 1 Laurent coefficients of the generating function, in the limit k — > oo it 
should be possible to extend (7) to more arbitrary functions, providing appropriate 
conditions for convergence are fulfilled. His last sentence reads "Only, before I knew 
what sort of conditions to impose on the generating function, we talked to Kakutani 
and Kakutani talked to Szego and the mathematicians got there first." 

Further explanation of that comment is given in [27], where he says that he 
had found "a general formula for the evaluation of Toeplitz matriceslj The only 
thing I did not know was how to fill out the holes in the mathematics and show 
the epsilons and the deltas and all of that" . Onsager adds that six years later the 
mathematician Hirschman told him that he could readily have completed his proof 
by using a theorem of Wiener's. 

Szego did publish his resulting general theorem [29], [30] on the large-size limit 
of a Toeplitz determinant, but not until 1952. 

So Onsager had a derivation of Ai in 1949, but looked for a more general way to 
calculate A^. For his rational generating functions he actually proved what is now 
Szego's theorem. He could see the extension to more general functions, and that it 
gave the same result (14. ip as his previous method. However, he lacked a rigorous 
proof of that extension and did not pursue the matter when the mathematicians 
became interested. 

The author also has a copy of a typescript, given to him by John Stephenson, 
entitled "Long-Range Order" and with three names hand-written at the top: 
"Onsager" nearest the title, "B. KAUFMAN" above that, and "R. B. Potts" above 
that. It shows how the magnetization can be obtained from eqn. 43 of [24], using 
Onsager's working of [25] and deriving the formula of the above-mentioned letter 
from Kaufman to Onsager. This must be a draft of a paper by Onsager and 
Kaufman on the subject. The author intends to make this typescript available 
at < http: / / arxiv.org/ > . 

| Refs. [26], [37] are reprinted in Onsager's collected works, pages 232 - 241 and 37 - 45, 
respectively. [35] 
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4-2. Yang and others 

The first derivation of Ai published was in 1952 by C.N. Yang, jl] He first 
wrote Ai in terms of the two maximal eigenvectors of the transfer matrix. He then 
used Kaufman's spinor operators to write Ai as the product of eigenvalues (i.e. a 
determinant) of an L by L matrix and evaluated the eigenvalues in the limit L — > oo. 

Later, combinatorial ways were found of writing the partition function of the 
Ising model on a finite lattice directly as a determinant or a pfaffian (the square 
root of an anti-symmetric determinant) . [3lT [32] . Then it was realised that the 
problem could be solved by first expressing it as one of filling a planar lattice 
with dimers.[33l E3 ESI |36] In 1963 Montroll, Potts and Ward[5] showed that the 
magnetization could be written as the ratio of two determinants, thereby obtaining 
a combinatorial proof of the result (14. ip . The numerator determinant is indeed 
Toeplitz and they evaluated its large-size limit by using Szego's theorem. 

5. Six- vertex model 

For any iV-state model (with nearest neighbour interactions) on a lattice of L 
columns, the transfer matrix T is of dimension N L . So is the hamiltonian if of a 
quantum mechanical system of iV-state spins on a line of L sites. Calculating the 
largest eigenvalue of T corresponds to calculating the lowest energy state of H. 

In 1966 C.N. Yang and CP. Yang [37] extended the Bethe ansatz by using it to 
calculate the ground-state energy of the hamiltonian of the anisotropic Heisenberg 
chain, also known as the XXZ chain. The following year Elliott Lieb used this 
method to calculate k for three particular models: ice, F and KDP.[6l [7J [8] These 
are all special cases of a more general zero-field "ice-type" or "six-vertex" model. 
The solution of this model was given in the same year by Sutherland. [9] Finally, the 
solution of the general six-vertex model was given by C. P. Yang [38] and Sutherland 
et al |39j. 

In these models one places arrows on the edges of the square lattice so that at 
each vertex there are two arrows in and two arrows out (this is known as the "ice 
rule"). There are six ways of doing this, as shown in the upper diagram in Figure 
[31 In general we assign to these six vertex configurations the weights u±, . . . ,uq, as 
in the figure. 

One can place a line on any horizontal edge bearing a left-pointing arrow and a 
line on any vertical edge bearing a down-pointing arrow; other edges one regards as 
empty. The six possible vertex configurations are then represented as in the lower 
diagram of Figure [31 using dots to represent empty edges. Separating the lines in 
the second vertex as indicated, it becomes apparent that these lines are continuous, 
and can be viewed as moving generally up and to the right through the lattice. 

If we impose cyclic boundary conditions, then it follows that the number of 
lines in one row of vertical edges must be the same as the number in the row above. 
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Figure 3. The six arrangements of arrows at a vertex in the six- vertex 
model . 



(Lines can move out at the right-hand boundary, but must re- appear at the left.) 
Hence if there are n lines in one row, then there are n lines in all rows. This is 
the "line-conservation" or "conservation of down arrows" property. The row-to-row 
transfer matrix T therefore breaks up into L+l diagonal blocks, one for each possible 
value (0, 1, . . . , L) of n. Two typical arrangements of lines in two adjacent rows are 
shown in Figure HI 



Vi 2/2 Vi 2/2 



Xi X 2 Xi x 2 

Figure 4. The two typical arrangements of lines in adjacent rows of the 
six- vertex model. The yi, . . . y n must interlace the x\, . . . ,x n . 
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5.1. The Bethe ansatz 

Here I shall briefly outline how the Bethe ansatz is applied to the usual row- 
to-row transfer matrix. In some ways it is simpler if one instead uses the transfer 
matrix in the diagonal SW to NE direction, since then lines only move at most one 
position at a time. Alternatively, helical boundary [ID] conditions mean that one 
only has to add a single vertex at a time. 

However, it is worth the effort to work with the row-to-row matrix, since the 
final equations simplify and can be compared with hamiltonian calculations such 
those mentioned above by Yang. Further, we shall find that they manifest the 
important commutation property of the transfer matrices. 

The method is discussed in jJTJ H2] and, for the zero-field case (when u\ = u 2 
and Cc?3 = W4), in sections 8.2 - 8.4 of |43j . Here I outline the extension to arbitrary 

U)\, . . . , U3q. 

Let A be an eigenvalue of the transfer matrix T and / the corresponding 
eigenvector. Each element of / is associated with a configuration of arrows or 
lines on a row of vertical edges of the lattice. Write f(x\, . . . ,x n ) for the element 
corresponding to lines (down arrows) at positions xi, . . . , x n be f{x\, . . . , x n ). Write 
the set xi, . . . x n as X, and y±, . . . y n as Y. We must have 

1 < xi < x 2 < ■ ■ ■ x L < L (5.1) 

and similarly for Y. 

Then the eigenvalue equation is 

A/PO = Y,C(X,Y)f(Y) + Y^D(X,Y)f(Y) , (5.2) 
where the first sum is over all yi, . . . , y n such that 

xi < yi < x 2 < y 2 < ■ ■ ■ < x n < y n < L , 
and the second sum is over 

1 < yi < x\ < y 2 < x 2 < . . . < y n < x n . 

Here X is the configuration of arrows in one row, and Y is the configuration in 
the row above, as in Figure [H The coefficients C(X,Y), D(X,Y) are the products 
of the Boltzmann weights of the vertices in between. Note that special cases occur 
when an X{ is equal to a yj, as then we can get a vertex of type 2 or type 3, with 
weight uj 2 or CJ3. We must also ensure that we do not include cases where two of the 
yi are equal, e.g. y± — x 2 — y 2 in the first sum. 

5.2. The case n = 1 

When n — 1 the equation ( 15. 2ft is (writing x±,yx simply as x, y) 

L 

Af(x) = uj^ l u 3 f(x) +u 5 u 6 J2 wf +x " y_1 w|" :E ~ 1 f(y) 

y=x+l 
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+ u 5 u 6 j2 u i~ v ~ luj 4 +v ~ x ~ 1 f(y) + uj i uj 4~ 1 f( x ) ■ ( 5 - 3 ) 

y=l 

We try the solution 

f(x) = z x 

and find the RHS of Q is 

(^C + u^M)z x + p(z)u x ut x (l ~z L ) , (5.4) 



where 



cji(o;i — u^z) 



uiu 3 + (u 5 uq-u 3 u 4 )z 

£ = £{z)= r 5.5 

Ui[U)i — u&z) 

M _ J\/l( z ) ^2 -uj 5 uj 6 - U 2 UjZ 

U^Ui — UJ^z) 

Note that M. is this section is not the same as the magnetization M. discussed 
elsewhere in this paper. It follows that (I5.3P is satisfied if 

A = u[C{z) + u±M{z) and z L = 1 . 

Thus there are L possible values for z (all lying on the unit circle), giving the L 
eigenvalues A in the n = 1 block of the transfer matrix. 

5.3. The case n = 2 



When n = 2 we first try 

f{x 1 ,x 2 ) = z Xl z*> . 

Then various terms come from the sums represented by the two diagrams in Figure 
HI Considered as functions of Xi,x 2 , the only ones proportional to zf 1 z 2 2 are 

~u[£{z{)C{z 2 ) + oo^M(z 1 )M(z 2 )] z x S z x * , 

so these terms in the eigenvalue equation are satisfied if 

A = uj[C( Zi )C{z 2 ) + u^M{z 1 )M{z 2 ) . 

We shall see that the contributions from the boundary terms lead to the 
"conservation of momentum" equation 

z[z L 2 = 1 . (5.6) 

Both these last two "energy" and "momentum" equations are unchanged by 
permuting z\ and z 2 , so we could more generally try an ansatz of the form 

f(x 1 ,x 2 )=A 12 z?z?+A 21 z?z? , 

where A± 2 and A 21 are coefficients that are at our disposal. Then the boundary 
terms contain factors p{z\) or p(z 2 ) and vanish if 



A 12 = z[A 2X , A 21 = z%A 



12 , 
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from which (I5.6P follows. 

The other terms that arise are "unwanted internal terms" terms from the first 
diagram of Figure H] containing a factor (ziz 2 ) X2 and from the second diagram 
containing a factor (ziZ 2 ) Xl . Both of these vanish if 

s(z 1 ,z 2 )A 12 + s(z 2 ,z 1 )A 21 = , (5.7) 

where s(zi,z 2 ) = UiU^Ai^z^C^) — oj 2 uj 3 . Removing factors that cancel out of 
(15. 7p . we can take 

S{Zi,Z 2 ) — . ^O.OJ 



5.4- Arbitrary n 

It turns out that these results generalise easily to arbitrary n, becoming 

f(X) = E^---C n . ( 5 - 9 ) 
p 

where the sum is over all n\ permutations P = {pi,p 2 , ■ ■ ■ ,Pn} of the integers 
{l,...,n} and Ap = A Plr „ )Pn . This is the Bethe ansatz. The unwanted internal 
terms give the equations 

s(z Pj , z Pj+1 )A P + s(z p . +1 ,z Pj )A P Q) , 1 < j < n (5.10) 

where ?W differs from P in the interchange of pj with Pj+i- There are many more 
equations than unknowns in (I5.10p . but they have the solution 

A P = e P J] s(z Pj ,z Pi ) , (5.11) 

l<i<j'<n 

where ep = +1 for even permutations, —1 for odd permutations. 
The boundary terms give 

Ap = z pi A P2 ,..., Pn , Pl (5-12) 
Eliminating all Ap between these two equations, we get the n "Bethe equations" 

zf = (-!)«-* ft • (5.i3) 

These imply z± ■ ■ ■ z^ = 1 and in general define Z\ , . . . , z^ . There are many 
solutions, corresponding to the different eigenvalues. Finally, the wanted terms 
in the eigenvector equation give the eigenvalue as 

A = u)fC(z 1 )---£(z n ) +u L i M{z l )---M{z n ) . (5.14) 

For the ice model, where uj\ = ■ • • = the full working is given in Ref. [6] 
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5.5. Free-fermion case 
Note that when 

OJ1OJ2 + UJ 3 U4 — OJ^UJq = 

the function s(zi,z 2 ) in (I5.8P is symmetric, so (15.13j) simplifies to 
4 = ("I)"" 1 , 

and the Bethe equations can be solved explicitly for all the eigenvalues of T. 

This is known as the "free-fermion" case. Like the Ising model, it can also be 
solved combinatorially by pfaflians. [H] 



5.6. Zero-field case 

The zero-field case is when 

0J\ = 0J2 = a , 0J3 = Ui = b , CJ5 = ojq = c , (5.15) 

a, b, c being positive parameters. The vertices 5 and 6 always occur in pairs, being 
sinks and sources of horizontal arrows, so there is no restriction in taking U5 = Uq. 
The first two conditions are definitely constraints. Then from (15. 8p . 

s(zm,Zj) = l-2Az j + z m z j , (5.16) 

where 

It turns out to be useful to express these equations in terms of new variables 
[i,v,(3j defined by 

a = sinv , b = sin(/i — v) , c — sin// , 

sin j3j 
sin(/i — (3j) 

This fixes the normalization of a, b, c and ensures that 



si n (y - /3 j - 11) M r z \ = s[n ( v ~ h + ^ 
sin(t> — (3j) sin(f — (3j) 



A = — cos// , C(zj) — — — — — , M(z 

From flBTTHjl . 

s{zm, Zj) = G mj sin(/3j - f3 m + fji) 

where G m j = Gj m = sin/i/[sin(/i — f3j) sin(/i — (3 m )]. The factors G m j cancel out of 
the RHS of (15.131) . so it depends on . . . f3 n only via their differences. 
If we define 

n 

(f)(v) = sin L v , q(v) = JJ sin(v - fa) , (5.18) 
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and note that A is a function of v, so write it as t(v), then f)5.13j) and (15. 14j) become 
(for j = l,...,n) 



where fjf = fi — ir. When — 1 < A < 1, we can choose fi,v to be positive, with 
< fi < n, < v < /I. 

For A < 1 and L even, the maximum eigenvalue lies in the block with 



and the Zi,...,z n lie on the unit circle, distributed about the point z — 1. In 
the thermodynamic limit of L large they form a continuous distribution and (15. 1 3[) 
becomes a linear integral equation for the distribution function. If we transform 
from the Zj to the (3j (the f3j lie on the vertical line Re(flj) = ji/2 in the complex 
plane), then this equation has a difference kernel and the integral is over a full period 
of the functions, so the equation can be solved by Fourier transforms, giving an exact 
expression for the free energy log ft. For the F model, Lieb found[7] that there was 
a transition of infinite order at A = —1, i.e. all derivatives of the free energy exist 
but there is an essential singularity of type e~ c ^ T ~ Tc ^ 2 , 7" being the temperature, 
T c the critical temperature and C a constant. For the KDP model he found[8] there 
was a first-order transition at A = 1 from a disordered phase (A < 1), to a frozen 
ordered state (A > 1) where the vertical (horizontal) arrows all point in the same 
direction. This was very different behaviour from that of the Ising model, which has 
a logarithmic singularity in the specific heat. 

5. 7. The model in a field 

The above working remains valid when electric fields are applied, making 
uj\ 7^ u>2 and/or u>s ^ LU4. The system is then either frozen or the maximum 
eigenvalue occurs when distributed along a curve in the complex plane. 

The resulting integral equation cannot in general be solved analytically by Fourier 
transforms. However, it can be solved numerically and at (and sometimes about) 
various special cases (including the free-fermion case) . [4"I] . [4"3| §8.12], [4"5 | |4T H W7 \ |4"8] 

5. 8. Matrix functional equations 

I had the good fortune to work with Elliott Lieb at the Massachusetts Institute 
of Technology from 1968 to 1970, and looked at a number of Bethe ansatz problems. 
In particular I considered the most general inhomogeneous six-vertex model that 
could be solved by the Bethe ansatz, using equations of the same form, but involving 
different parameters at different sites of the lattice. [49] In particular, instead of the 
simple exponentials z x - above, one uses a more general "single-particle" function 



t(v)q(v) = 4>{v) q(v — //) + 4>{jJL — v) q(v + jj!) 



(5.19) 
(5.20) 



n = L/2 
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In 1970 my wife Elizabeth and I left MIT to return to Australia via England, 
where we would spent time with my parents in Essex. We intended to make the 
journey from England to Australia by ship. Already it was unusual to do this rather 
than travel by air (a change that had occurred only in the previous decade). We 
had the choice of two ships, giving us either two or five months' break in England. 

We chose five months, and it was a good decision. Towards the end of that 
time I looked again at my MIT paper and realised that there were parameters in 
the final definition of the model that entered the eigenvalues of T, but not the 
eigenvectors. Indeed, this was apparent from Lieb and Sutherland's work on the 
zero- field six- vertex model: from (15. 9p - (I5.17p . the eigenvectors / depend only on 
the single parameter A = — cos \i. They are the same eigenvectors as those found 
by Yang and Yang [37] for the hamiltonian of the anisotropic Heisenberg chain. 

The Boltzmann weights of the model depend also on v, If we write the transfer 
matrix T as T(v) and keep A and fi fixed, then under quite general conditions it 
follows that two transfer matrices T(u), T(v), with the same value of A, commute: 
T{u)T{v) = T{y)T{u) . (5.21) 

There is a common eigenvector matrix V, independent of v, such that V~ l T(v)V is 
diagonal, its entries being the eigenvalues t(v). 

Form a matrix Q(v) such that V~ 1 Q(v)V is also diagonal, with entries the 
corresponding q(v) of (I5.20p . Then T(v) and Q{v) must also satisfy (I5.20p . i.e. 

T(v)Q(v) = <f)(v) Q(v - fjf) + (f}(n - v) Q(v + fjf) , (5.22) 

the <fi factors being scalars; Q(v) must commute with T(u), for all u, v. 

This matrix relation, together with some elementary observations, is equivalent 
to the Bethe ansatz. Together with the commutation properties, it implies (I5.20p . 
Since the elements of T are sums of products of L Boltzmann weights, each equal 
to a, b or c, and vertices 5 and 6 always occur in pairs, each element must be of the 
form e~ iLv times a polynomial of degree L in e 2lv . If the elements of Q have a similar 
structure, but with L replaced by n, then the same must be true of the eigenvalues 
t(v),q(v); and q(v) must have n zeros (3i,...,(3 n . Setting v = (3j in (I5.20p (and 
noting that t(v) is entire, so finite), we obtain the Bethe equations (I5.13p . 



6. Eight-vertex model 

This was a new way of solving the six-vertex model. I knew from conversations 
with Elliott Lieb that an interesting generalization of the six-vertex model was the 
eight-vertex model, where one allows two more vertices, one with all arrows in, the 
other with all arrows out, with weights Uy,Ug, respectively^ 

§ Again, there is a special "free-fermion" case, when + — — ^7^8 = 0, which can 
be solved by pfafhans.fjl] 
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This general eight-vertex model does not have the "conservation of down 
arrows" property that was necessary for the Bethe ansatz, but an obvious question 
was whether it has the commuting transfer matrix property. If so, then could one 
extend the above T, Q relation (I5.22p to include this model, and so calculate k and 
the free energy? 

It turned out that the answer to both questions was yes, provided we restrict 
consideration to the zero-field model, with 

uj\ = 0J2 = a , Cc>3 = Co>4 = b , = ujq = c , uj-j = us = d . (6.1) 

(The last two equalities are not restrictions, as uj-j and u$ must occur in pairs, as do 
and uq.) If we define 

2{ab + cd) ab + cd 

then the transfer matrices of two models with different values of a, b, c, d, but the 
same values of A and T, commute. 

Sutherland had shown in 1970 [50] that each transfer matrix commutes with the 
hamiltonian of the XYZ chain (with coefficients A,T), which implies the required 
commutation relation between transfer matrices. I was unaware of this, but verified 
it myself and went on to construct the matrix Q(u). This also satisfies (I5.22p . but 
now 4>{u) is a product of Jacobi elliptic theta functions, as must be the eigenvalues 
T(u), Q{u). They are doubly periodic in the complex w-plane and entire; 0(ii) 
and T(u) have L zeros per period parallelogram, while Q(u) has L/2. Again, 
f)5.22p is sufficient to determine the eigenvalues T(u),Q(u). The results for the 
critical exponents were unusual and excited interest - they vary continuously with 
a parameter /1 of the model that corresponds to the fi of the six-vertex model. For 
example, the critical exponent a of the specific heat is 

71 

a = 2 - - . 
/1 

The eight- vertex model includes the Ising and six- vertex models as special cases. 
The Ising case is /1 = 7r/2, A = 0, corresponding to the logarithmic singularity; 
/i = 0, A = —1 corresponds to the F model transition, with a = —00; and fi = tt, 
A = 1 to the KDP model transition, with a = 1. 

Of course, all this took many weeks to work out, and some of it was done on 
the ship from England to Australia. I shall always wonder if I would have had this 
idea if I had not had that five months' break in England - it can be a good idea to 
give the mind time to relax and go on auto-pilot. 

As with the Ising model, it was harder matter to obtain the spontaneous 
magnetization and spontaneous polarization. I obtained the polarization in 1973 for 
the six- vertex model in the anti-ferroelectric regime A < —1 by a direct calculation 
of scalar products of Bethe eigenvectors: the method was somewhat tortuous, as 
is indicated by the fact that the proof of one step depended on taking n to be a 
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prime number! Michael Barber and I conjectured the spontaneous magnetization of 
the eight- vertex model in that year, [51] and Stewart Kelland and I the spontaneous 
polarization in 1974, [52] but derivations had to wait for the development of the 
corner transfer matrix method in 1976. [53] 

7. Yang-Baxter relation 

Two-dimensional nearest-neighbour lattice models can be grouped in three 
classes: 

• Spins live on sites and interact along edges (e.g the Ising model) . 

• Spins live on edges and interact at a vertex (e.g. the six and eight-vertex 
models). 

• Spins live on sites and all four spins (on the square lattice) round a face interact 
(the eight- vertex model can also be formulated this way). 

The condition for two transfer matrices to commute is a local one, involving 
the Boltzmann weights W of particular edges, vertices or faces. It takes different 
forms[I7], (121 eqns. 9.68, 13.3.1] for the three above cases, and these can be 
represented graphically as in Figure [5j In each case W\ can be interpreted as the 
weight (edge, vertex or face) of one transfer matrix, Wi of the other, and as a 
supplementary weight in the equation. Each figure represents the partition function 
of a small graph. The outer spins (on the open circle sites or the six exterior edges 
of the second figure) are fixed and the inner ones (solid circles or the three edges of 
the inner triangles) are to be summed. The partition functions on each side must 
be equal, for all values of the corresponding exterior spins. 

The first form is the star-triangle relation, which was used by Onsager.pQ The 
second was used by McGuire.[5l] 

7.1. Other models 

This general technique turned out to be very useful in identifying and solving 
statistical mechanical models: some two-dimensional examples are 

• Hard-hexagon model in 1980. [H] 

• Fateev-Zamodchikov model in 1982.[i~2] 

• Kashiwara-Miwa model in 1986. [13] 

• Solvable Y-state chiral Potts model in 1988. [IT] 

In addition, the Yang-Baxter equations can be straightforwardly extended to 
three dimensions. The first solution to the resulting "tetrahedron" equations was 
given by Zamolodchikov in 1981. [H] Further solutions were found later, e.g by 
Bazhanov and Baxter. [15J As a general rule these three dimensional models are 
critical and their weights are not necessarily real and positive. This contrasts with 
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Figure 5. The three forms of the Yang-Baxter relation. 

(say) the hard hexagon model, which is a good model of a triangular lattice gas, 
exhibiting the different phases, and with critical exponents that can (because of 
universality) be compared directly with experiment. [55] 

7.2. Mathematical difficulties 

Of all the models mentioned, the Ising model is undoubtedly the simplest. 
Because of its underlying free-fermion structure, all the eigenvalues can be evaluated 
explicitly, even for a finite lattice width L. This contrasts with the other models, 
where one usually does not know any explicit solutions of the non-linear T, Q (or 
equivalent) relation, and the best one can do is to use analytic techniques to calculate 
the largest and near-largest eigenvalues of T in the limit of large L. In general there 
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is no simple direct product property for the eigenvalues. 
7.3. Rapidities 

For all these two-dimensional solvable models, the Boltzmann weights W (and 
W) depend on both the spin variables (e.g. a, b) and on other "rapidity" variables 
p,q: 

W = W{a,b) = Wpq(a,b) , W = W(a, b) = W pq (a, b) . 

These rapidity variables are associated with the dotted lines in Figure [TJ (For the 
moment ignore the break in the line below a\.) Each edge of the original lattice of 
solid lines is intersected by two dotted lines, one horizontal and one vertical. The 
vertical (horizontal) dotted lines have rapidities p or p' (q or q'). In general each line 
may have its own rapidity variable, but for an homogeneous model all the vertical 
rapidities p must be the same, and similarly for the horizontal rapidities q. 

The rapidity variables can be chosen so that two row-to-row transfer matrices, 
with different values q, q' of q, commute: 

T q Tqi = T q i T q . 

If the weight W\ in Figure [5] depends on two rapidities q and r, and W% on r and p, 
then W3 depends only on p and q: W3 is independent of r. 

The effect of the Yang-Baxter relation is that one can move these rapidity lines 
around without changing the partition function. [56] If these moves do not cross the 
spin o\ in Figured! then the RHS of (12.31) is also unchanged. If the lattice is infinitely 
big and o~\ deep within it, then any particular rapidity line p (or q) can be moved 
infinitely far from a%, where we do not expect p (or q) to contribute to the RHS of 
(14. ip . It follows that the magnetization M. must be independent of all the rapidity 
variables, and indeed that is implied by the result (14. ip . (For the Ising model, K 
and K depend on the rapidities, but k is a constant, the same for all edges of the 
lattice.) 

7-4- Rapidity difference property 

With the exception of the solvable chiral Potts model, all the two-dimensional 
models mentioned above have the rapidity difference property, i.e. their weights 
W, W depend on the corresponding rapidity variables p, q only via the difference 
p — q (the variable u in (15.221) is such a difference). This is significant. It means 
that W, W are trigonometric or elliptic functions of p — q. In the limit of a large 
system the partition function per site k can be obtained easily by the "inversion 
relation" trick [57] and the spontaneous magnetization by the corner transfer matrix 
method. [I3j chapter 13] These methods do involve assumptions about the analytic 
properties of the variables, for instance that k is an analytic function of u = p — q 
in some vertical strip in the complex w-plane. 
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8. Solvable Chiral Potts model 

For the solvable chiral Potts model, W(a, b) = W(a—b) and W(a, b) = W(a— b), 
where 

fh\ n ^^P W[n) = { )n jj - <* x q 

Here 

u = e 2m / N , (8.2) 
fi p , x p , y p are three complex variables related by the two equations 

kx N p = 1 - h>/$ , ky? = 1 - k'np , (8.3) 
and k, k' are constants satisfying 

k 2 + k' 2 = l . (8.4) 

We can think of p = x p , y p } as a point on a three-dimensonal curve. This is 
the rapidity p. Similarly for q. If x p , x q , y p , y q lie on the unit circle and are ordered 
anti-clockwise in that sequence, then we can choose fi p , fi q so that the Boltzmann 
weights W, W are real and positive. The conditions (I8.3P for p, q ensure that 

W(n + N) = W(n) , W(n + N) = W(n) . (8.5) 

The model is therefore Z^-symmetric. It is chiral because W(—n) ^ W(n) and 
W(—n) ^ W(n), so reflection-symmetry (left-hand equals right hand) is broken. 

When N = 2 this is the Ising model, with e~ 2K = W(l), e" 2l? = W{1) 
(normalising the weights in (13. ip so that W(0) = W(0) = 1). 

For N > 2 this model does not have the "rapidity difference" property. 
This makes it much harder to calculate the thermodynamic properties. Even so, 
differential equations were written down [58] in 1988 that in principle define k and 
the critical exponent a was thereby found to be 1 — 2/N. An explicit expression as 
a double integral was obtained in 1990. [591 EOl EI] 

It is natural to define the magnetization as 

M r = (u r ^) = Z- 1 ^^ 1 n^,^) (8-6) 

for r = 1,. . ., N - 1. 

In 1989 Albertini et al conjectured [62] that for < k < 1 

Mr = {1 _ k 2y(N-r)/2N* ? (g _ 7) 

which is a beautifully simple formula that fits the N = 2 Ising case and the available 
series expansions. For k > 1 the system is disordered and the magnetization 
vanishes: Ai r = 0. 

It was very difficult to prove this conjecture. In his collected works C.N. Yang 
says his 1952 calculation of the Ising model magnetization took six months and was 
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the longest in his career. [55] The chiral Potts magnetization took much longer: it 
was not till 2005, after spending much time on and off mulling over the problem, 
that I succeeded[6l] in deriving the formula f!8.Tj) . 

8.1. Broken rapidity line derivation of Ai r 

The technique I used was invented by Jimbo, Miwa and Nakayashiki. |65j One 
breaks one of the dotted rapidity lines adjacent to the central spin <j\ and gives the 
two halves different values q\ and q 2 of the rapidity variable, as in Figure [TJ 

The effect of this is that one cannot remove these two half-lines away from o~\. 
One can still remove all the other rapidity lines, so now, using the definition (18. 6p , 
Ai r = function only of k, q± and q2 ■ (8.8) 

This is a generalization of the magnetization. 

However, one can rotate the two half-lines and cross them over. This gives 
functional relations for the generalized Ai r - For the models with the rapidity 
difference property, Ai r can only depend on qi, q 2 via their difference. If one makes 
plausible assumptions about the analyticity properties of this function (e.g. analytic 
in a particular vertical strip), then one can solve the functional relations and obtain 
Ai r - This provides an alternative to the corner transfer matrix method of calculating 
single-spin correlations. 

Again, however, life is more difficult for the chiral Potts model. It does not 
have the difference property and (ignoring the dependence on k, which we regard 
as a constant) Ai r is a function of two variables, not one. It is not obvious how to 
solve the functional relations, and what additional information is required. 

What I observed in 2005 was that if one took 

Xq 2 = Xq ± , ]Jq 2 = Uiy qi , fiq 2 = flq 1 , ($-9) 

and again made a plausible analyticity assumption, then one could evaluate A4 r for 
this case. 

Of course one wants it rather for the case when q 2 — qi, but since it is then a 
constant, independent of qi, it is sufficient to obtain it at the intersection of these 
two cases, when 

y q2 = y qi = . (8.10) 

In this way I was able to verify the conjecture ( 18.71) . I still do not know what the 
generalized Ai r is for arbitrary q±,q2- 

9. Superintegrable Chiral Potts model 

We return to the usual situation, where there are no broken rapidity lines and 
the magnetization Ai r of (18. 6p is independent of the rapidities. 
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There is a special "superintegrable" case of the chiral Potts model, when the 
vertical rapidities alternate, as in Figure^} taking the values p,p',p,p', . . ., where 

ay = y P , y P ' = x P , Hp' = iM> • (9- 1 ) 

Since Ai is independent of the rapidities, its value for the superintegrable case is 
the same as that of the general solvable chiral Potts model. 

The chiral Potts model may be the most difficult of the two-dimensional solvable 
models, but its superintegrable case has some remarkable simplifications. In fact if 
we impose cylindrical boundary conditions as in Figure [TJ with the top and bottom 
spins fixed to zero, then it has properties similar to those of the Ising model. (For 
N = 2, the superintegrable case, like the general solvable the general zero- 

field Ising model.) 

With these boundary conditions, the partition function Z is given by (12. 5p and 
M r by 

u^T n S r T n 'u . . 

M r = j , (9.2) 

where n is the number of rows below o\ and n' the number above, so n + n' = M. 
The vector u is defined by f)2.6p . Also, S r is a diagonal matrix with entries 

= c^n^Oi) ( 9 - 3 ) 

i=l 

again writing the spin set {01, . . . gl\ as s. 

9.1. Partition function Z 

Define a set of vectors u , . . . , itjv-i with elements 

K) s = f[S(o- u a) . (9.4) 

8=1 

Then u = uq, where u is the vector above, defined by (12.61) . 
Let R be the spin-shift matrix with elements 

L 

Rs,s' = n%i.^+ i ) • ( 9 - s ) 

i=l 

Then R N = 1 and the eigenvalues of R are l,ou, . . . ,a; Ar ~ 1 . Let Vq be the vector 
space of all vectors v such that Rv = u®v. Then u a = R a u and 

N-l 

v Q = N- 1 ' 2 u~ Qa u a E V Q (9.6) 

a=0 



and 



where 



Z = Zo + ^ + .-. + Zjv-i , (9.7) 
Z Q = Z Q {M) = VqT m vq . (9.8) 
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These matrices are of dimensions N L and Vq is of dimension iV L_1 . However, 
if we repeatedly multiply Vq by T, we generate a smaller sub-space, [HI [66] of 
dimension 2 m , where 



rn 



m Q 



(N-1)L-Q 

N 



(9.9) 



and we write [x] for the integer part of x. Further, we can choose the basis vectors 
(independently of q and k) so that 

■■®t m (9.10) 



and 



/ 



V Q 



A 



i 



v / 



A 



(9.11) 



V / 



Here each tj is a two-by-two matrix. So to calculate Zq we only need T in a 2' m - 
dimensional sub-space, within which it is a direct product, as in the Ising model. 
It follows that Zq = Zq(M) is a product of m simple factors. 



9.2. Magnetization 
From fl£2} and 
M r = 



JV-l,r 



(9.12) 



Zo + ■ ■ ■ + Zft/_1 

where 

W Qir = v Q T n S r T n 'v Q+r . (9.13) 

When JV = 2 we regain the Ising model and we know then[U |5l [19] for 
finite L that Wq^ t can be written as a determinant. We have just seen that the 
superintegrable chiral Potts resembles the Ising model in that Zq is a simple product. 
We therefore ask whether our Wq^ can be written as a determinant for all Nl 

The answer is yes. It can be written as a determinant. The working is given in 
a series of papers (not necessarily in logical order) . [20j ED, [22} [23], [70] 

Setting m! = mQ +r , we find that 

W Qtr = Z Q (n)Z Q+r (n')det[l + AA T } , (9.14) 

where A is an m by m! matrix. 

This is huge progress. If n = 3 and L = 15, then m = m! = 9. We have reduced 
the problem from one calculating the elements of powers of the transfer matrix T, 
of dimension 14,349,907, to one of evaluating a nine by nine determinant! 

Even so, we do want to take the large lattice limit. It is easy to allow n,n' to 
become infinite - we merely take the limit of the elements of the matrix A. We then 
want to let L — > oo. 
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In this limit even m and m' will become infinite, so we need a way of evaluating 
the determinant in (I9.14p . 

This was a tricky problem and I mulled over it for more than a year. If m = m', 
then A is square and invertible, so 

det[l + AA T ] = det A det[(A T )" 1 + A] . (9.15) 

The elements of A are of the form 

A l3 = , (9.16) 

where known parameters. 

The matrix A is therefore Cauchy-like and its determinant is [6 



, , a _ XX 'Ul<i<j<m(Ci - Cj)(Cj - c-) 

n^n™^-^) ' (9 ' 17) 

where X = Y[ i Xi 1 X' = ]Jj x'y Further, the inverse of the transpose of a Cauchy-like 
matrix is also Cauchy-like, the elements also being of the form (I9.16p . with only 
Xi,x'j changed to some values y^y'y 

The elements of the desired sum are therefore 

UA T y l + A) = VlV ' 3 + X fi . (9.18) 

v ' l i ci — c\j 

This is not a Cauchy-like matrix, rather it is akin to a Pick matrix of 
displacement rank 2. [68] There are fast computational algorithms for numerically 
calculating the determinants of such matrices. [691 Chapter 1] but I know of no 
explicit expression for the answer. 

However, I finally looked at the numerator in the desired limit n, n' — > oo. It is 
a rational function of Aj, A^, where 

d = (1 + k' 2 - Xi 2 )/2k' , c' 3 = (1 + k' 2 - X'/)/2k' (9.19) 



so 

Cl - c' 3 = (Af - A, 2 )/2^ . (9.20) 
Further, the numerator of (19. 18ft turns out to be of the form 

yiy'j + x i x' j = Sis'^Xi + X'j) (9.21) 
so the factor Aj + A^- cancels out of (19.181) . leaving 

((AT) - 1 + 4 ^^M ( , 22) 

This is a Cauchy-like matrix and one can obtain its determinant from the general 
formula (19. 17[) . In this way one can obtain WQ, r /ZQ{n)ZQ +r {n') as a double product 
over i,j — l,...m, for finite L. 

In general m, m' differ by at most one. If they are different, one can add a row 
or column to A so as to make it square and still Cauchy-like, while leaving 1 + AA T 



Developments in exact solutions in statistical mechanics 



25 



unchanged. The same cancellation then occurs in (19.181) and again one obtains the 
double product expression for Wq, t . 

Finally one lets m, m' — > oo and evaluates the double product, using a Wiener- 
Hopf factorization. We again of course obtain the formula (18.71) . 

So this provides a rigorous proof of the magnetization of the superintegrable 
chiral Potts model, and hence of the Ising model. It does not involve Szego's 
theorem. The full calculation is given in refs. [19] - [23] and [70] . A related derivation 
of the magnetization via the pair correlation function is given in [TTj . 

10. Summary 

There have been many developments in the statistical mechanics of lattice 
models since Onsager's famous solution of the Ising model in 1944. The general 
trend has been one of increasing complexity, first to models without the direct 
product transfer matrix property, then to the chiral Potts model where one loses 
also the useful rapidity-difference property. 

But if one uses cylindrical boundary conditions with fixed spins on the top and 
bottom rows, then the wheel comes full circle with the superintegrable case of the 
chiral Potts model. One obtains determinantal expressions for the magnetization, 
analogous to those found for the Ising model by Onsager|24j. Yang[4] and Montroll, 
Potts and Ward [5]. 

One can evaluate these determinants, not by calculating matrix eigenvalues, 
nor by using Szego's theorem, but by identifying the determinants as Cauchy-like. 
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